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Abstract 

Small- World Networks (SWNs) represent a fundamental model for the comprehension of many complex 
man-made and biological networks. In the central nervous system, SWN models have been shown to 
fit well both anatomical and functional maps at macroscopic level. However the functional microscopic 
level, where the nodes of a network are composed by single neurons rather than brain areas, is still 
poorly understood. At this level, although recent evidences suggest that functional connectivity maps 
exhibit small-world organization, it is not known whether and how these maps, potentially distributed 
in multiple brain regions, change across different conditions, such as spontaneous activity and effective 
sensory stimulation. 

We addressed these questions by simultaneous multi-array extracellular recordings in three brain regions 
diversely involved in somatosensory information processing: the ventropostero-lateral thalamic nuclei 
(VPL), the primary somatosensory cortex (SI) and the centre- median thalamic nuclei (GM). From both 
spike and Local Field Potential (LFP) recordings, we estimated the hmctional connectivity maps by 
using the Normalized Gompression Similarity (spikes) and the Phase Synchrony (LFPs). Then, by using 
graph-theoretical statistics (transitivity, characteristic path length, small- worldness), we characterized 
the functional maps topology both diiring spontaneous activity and sensory stimulation. 
Our main results show that: (i) A SWN organization is present in spikes and LFPs during spontaneous 
activity; (ii) After stimulation onsets, while substantial functional map reconfigurations occur both in 
spike and LFPs, small-worldness is nonetheless preserved (iii) At the LFP level the stimulus triggers 
a significant increase of inter-area functional connections without modifying the topology of intra-area 
connections; (iv) In computational simulations on large-scale networks, the values of small-worldness 
statistics are consistent with those obtained in experiments. 

Our results suggest that coherent activity among neurons from multiple areas allows for the integration 
of local computations occurred in distributed functional cell assemblies according to the principles of 
SWNs. 

Author Summary 

Neuronal assemblies, sequences of neuronal activations, seem to represent the functional unit of informa- 
tion processing. However, it remains unclear how groups of neurons may organize their activity, during 
information processing, working as a sole functional unit. One prominent principle in complex network 
theory is covered by small- world networks, in which nodes are easily reachable by each other but preserv- 
ing low connection densities. Small-world networks are already found in networks observable in human 
or primate brain areas while their presence at the neuronal level remains widely debated and unclear. 
The aim of this work was to investigate the possibility that neuronal assemblies, encompassing multiple 
brain regions, can work as small-world networks. We investigated the coherent neuronal activity among 
multiple rat brain regions involved in somatosensory information processing. We found that the recorded 
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neuronal populations represented small-world networks and that this organizing principle did not change 
during stimulations. These results were also confirmed on large scale computational models. 
This work suggest that coherent activity of neurons from multiple brain areas, which can compose in- 
dividual functional units, promotes the integration of local computations, the functional principles of 
small- world networks. 



Introduction 

Neurons in the brain form highly composite networks sustained by a complex and variously distributed 
thread of synapses. Although the detailed anatomical connections are still under investigation |lj it 
has been shown that, during in-vivo recordings, active neurons form functional assemblies that do not 
necessarily depend on the underlying anatomical connectivity [2]. Therefore, while anatomical connec- 
tivity is stable for relatively long times, the functional connectivity is highly dynamic and depends on 
the particular tasks triggered by internal and external events. Critically, the organization principles that 
control the functional connectivity among single neurons and small neuronal assemblies are still poorly 
understood [3j. The early hypotheses |4 that distributed groups of neurons can operate as a functional 
unit through coordinated activities have now been confirmed in observations of transient cell assemblies 
over different cortical regions [5|. It remains however unclear how these dispersed neurons may gather in 
a functional unit for information processing. 

In the last 15 years, thanks to a substantial advancement in the complex network theory. Small- World 
Networks (SWNs) emerged as a paradigmatic model and provided the analytical tools to explain a large 
set of complex networks in the most diverse scientific areas [6j[7]. Furthermore, it has been proved that 
small-world networks represent optimized structures accomplishing adequate balance between informa- 
tion transfer efficiency and reliability [7-11 . Following this trend, large networks of brain areas, when 



studied from a macroscopic perspective by imaging techniques, have been shown to functionally organize 



as SWNs 12 13 . Contrarily, the functional connectivity emerging at the microscopic level of single 



neuron networks is still poorly understood, although recent works provided evidence of SWNs in small 



local neuronal populations 14-16 . Specifically it is not clear whether functional groups of neurons from 
multiple brain regions display a SWN organization. 

In the present work we address these points by means of spike and LFP simultaneous multi-array record- 
ings from the ventropostero-lateral thalamic nucleus (VPL), the centro-median thalamic nucleus (CM) 
and, the primary somatosensory cortex (SI) [T7||l9] . Our results confirm the presence of SWNs during 
spontaneous activity and provide evidence that, during stimuli, small-world organization principles are 
preserved in spite of massive topological reconfigurations. 



Consistently with other previous findings 5p0 ■ 22 , our results provide evidences for a model of functional 



organization where distinct functional neuronal assemblies are sparse and encompass multiple brain areas. 



Results 

Functional Connectivity in Spiking Activity 

Our first aim was to estimate the topology of functional connectivity graphs obtained from simultaneous 
recordings of spiking activity from neuronal populations in VPL, CM and SI. Because specific tactile 
stimulation can potentially exert a dramatic effect on functional connectivity, we set out to evaluate 
functional connectivity graphs both during spontaneous and evoked activity. In order to compute the 
functional connectivity graphs in these two conditions, we performed pairs of 10 minutes recordings, 
the first without any stimulus, the second by delivering fast tactile pulses (see methods) at randomized 
intervals (500 ± 200 ms) on the five different digits. We estimated the functional connectivity in neuronal 
spiking activity by using the Normalized Compression Similarity (NC'S) function to detect functional 
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couplings between pair of neurons. NCS is defined in the range [0, 1], where indicates no interaction 
and 1 indicates an exact correspondence between the firing patterns of the two neurons considered. This 
measure was selected instead of more classical ones because of its powerful property of capturing both 
short-range (synchronous) and long-range (delayed) interactions between neurons (see Methods, Fig. 
E). Thanks to this NCS property we could scan along the recorded activity in search for functional 
relations by using sliding windows of different durations (50, 250, 500, 1000 ms). The width of a slid- 
ing window defined the largest possible delay at which an interaction could be detected (see Methods). 
After the NCS estimation on all possible pairs of simultaneously recorded neurons, the resulting ad- 
jacency matrices were thresholdcd and binarizcd in order to obtain the functional connectivity graphs. 
Finally the small- world statistics C (transitivity), L (characteristic path length), S (small- worldness) 
were estimated on these graphs. The first measures the tendency of neurons to segregate into separate 
sub-networks, the second measures the average path length between nodes and the third is computed as 
S = [C /C^) / [L/ V). The terms C and U' represent normalizing values estimated on Erdos-Renyi ran- 
dom networks [23| with the same number of connections observed in the original graphs. A small-world 
network emerges whenever S > 1. 

First, we measured the small- worldness statistics S on our spontaneous activity recordings and we found 
that it was significantly larger than unity (P < 0.002, nonparametric Wilcoxon ranksum test), irrespec- 
tive of the time window used for detecting the functional connections (50, 250, 500, 1000 ms; Fig. [2]). 
Thus S satisfied the small-word conditions with invariance on different time timescales and the average 
values of C and L were respectivly 3.53 ± 0.4 (SD) and 0.57 ± 0.06 (SD). Larger time windows reported 
no small- world network configurations. 

Then we compared the average functional connectivity in time windows of 100 ms duration before and 
after the stimulus onset. The functional connectivity graphs underwent significant topological reconfigu- 
rations after the tactile stimulus delivery. As an example, in Fig. [3] the neuron 32 switches from a fully 
disconnected state (Fig. |3]^, relative to pre-stimulus condition) to highly connected (hub) state (Fig. 
[3J3) after the stimulus onset. Surprisingly, notwithstanding these profound changes, all the small-world 
statistics computed on the pre- and post-stimulus functional connectivity graphs were not significantly 
different. This can be seen in Figure |4 were we compared, along the five digits, C, L and S before 
(respectively [4|\, [4]C, |4|i;) and after ([4|3, 4|D,|4j^') the stimulus onset. Therefore the average small- world 
properties were not perturbed by tactile stimulations. 

Furthemore, we performed a correlation analysis between the quantitative changes in functional connec- 
tivity and the firing rate variations induced by tactile stimuli. We found that the evoked changes in 
functional connectivity graphs were significantly greater (P < 0.005, ranksum test) than in spontaneous 
activity. This means that stimuli effectively modulated the functional connections among neurons. In 
fact, the firing rate of neurons were positively correlated with the evoked connectivity changes {R = 0.637, 
P < 0.001, t-test; Fig. SI). Namely, whenever the recorded neurons were effectively involved in the re- 
sponse to stimuli, these stimuli recruited concurrent functional connection changes proportional to the 
evoked firing rates. 



Functional Connectivity in LFPs 

Our second aim was to estimate the topology of functional connectivity graphs obtained from the LFP 
activity recorded simultaneously to the spiking activity. LFP activity mainly refiects pre-synaptic events 
and therefore its investigation provides novel information that integrates the results on (post-synaptic) 
spiking activity reported in the previous section. We estimated the LFP functional connectivity maps by 
following the same sequence of analyses used for spiking activity. First, we estimated the functional con- 
nectivity between all possible channel pairs in order to generate the adjacency matrix, then we binarizcd 
this matrix by using a threshold (see Methods) and finally we computed the small-world statistics C, L 
and S. 

In order to estimate the functional connectivity between pairs of channels we used the phase synchrony 
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measure 7 (see Methods) in place of NCS. Phase synchrony is more suitable than NCS for continuous 



signals and it has been widely applied to EEG and LFP analyses 24 . It is normalized in the range [0, 1], 
where unity occurs when phase coupling is exact and zero when it is null. 

On spontaneous activity recordings, LFP functional connectivity graphs exhibited time-invariant small- 
world properties as the S statistics was consistently larger than unity (Fig. [5]A_-D). The average values of 
C and L over the different time-scales were respectively 2.98 ± 0.3 SD and 0.64 ± 0.05 SD. As a control, 
we computed C and U" by generating Erdos-Renyi graphs with the same number of connections and we 
found that their values were significantly smaller (P < 0.0005 for both comparisons, ranksum test). 
On evoked activity recordings, stimulus occurrence triggered significant topological reconfigurations in 
functional graphs, as shown, as a representative case, in Fig. [6] The picture reports a functional con- 
nectivity graph estimated before (Fig. |6]^) and after (Fig. [6j3) the stimulus onset. During spontaneous 
activity, CM (green), VPL (blue) and SI (red) showed tight intra-arca synchronization and poor inter-area 
synchronization. After the stimulus onset inter-area synchronization emerged while intra-area synchro- 
nization was maintained roughly constant. 

We asked whether these observations on a single recording could generalize to our full dataset. In order 
to test for this possibility we first divided our LFP recordings into two classes: responsive and non- 
responsive. LFP recordings were considered responsive when the average evoked firing rate, measured 
from the same recording channels (see Methods), was larger than the mean -I- 5 SD of the basal firing 
rate. We found that the number of intra-area (local) functional connections was preserved across the 
three conditions of spontaneous, non-responsive and responsive LFP activity (P = 0.63 for CM, VPL 
and SI, ranksum test). However, during responsive LFP recordings, the average number of inter-area (or 
global) functional connections was substantially larger for all possible inter-area combinations (P < 0.002 
for CM- VPL, CM-SI, VPL-Sl, ranksum test). 

Finally we compared pre- and post-stimulus functional graphs by computing the C, L and S statistics 
in 50 ms windows before and after the onset of tactile stimuli. We found that both conditions exhibited 
small- world properties and the statistics were not different (Fig. [?]). Interestingly, the tactile stimulus 
onset triggers a substantial increase in the number of inter-area connection but does not have a significant 
effect on the intra-area connections (Fig. S2). 

We further investigated the possibility, by NCS function, that asynchronous LFP phase patterns can 
however support small- world functional configurations. We found that no small- world networks emerged 
from such analysis (Fig. S3) suggesting that LFPs require tighter time constraints than spiking activity 
to put in place information processing. 

We can conclude that small-world statistics are preserved in LFPs during both spontaneous and evoked 
activity. 

Functional Connectivity in Simulated Spiking Activity 

In the previous sections we provided evidences, both at pre- and post-synaptic level, that cortical and 
subcortical neurons are functionally organized as small -world networks. Our results are consistent with 
recent findings from large scale cortical recordings 9 T0p!2 and suggest the existence of sparse functional 



networks composed of neuron assemblies that encompass multiple cortical and subcortical areas; impor- 



tantly, each functional network would recruit its units not necessarily in close proximity [5j|20 21 . 
Building on these results we conjecture that the range of values we estimated for S in our electrophysio- 
logical recordings critically depends on the number of underlying functional networks that are active at 
the time of recordings. This number represents a hidden variable and cannot be accessed experimentally. 
Therefore, in order to understand the dependence of S on the number of functional networks, we relied 
on a simulation approach. 

First we implemented a large-scale simulation of N neurons that embeds k functional networks. Then 
we estimated S by systematically varying k. We simulated a neural network of iV = 100000 neurons, 
about five fold the estimated size for a rat cortical column in the somatosensory system (25l. The size 
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of each functional network was chosen by samphng from a discrete uniform distribution ^ U{50, iV/100). 



This distribution and its parameters are consistent with the results of 26 . For further details about the 
implementation see the relative Methods section. 

At each run we randomly extracted 100 units (comparable with the number of units that we could si- 
multaneously record during the experiments) and we computed the S statistics. We observed that S 
underwent a bimodal behaviour (Fig. |8j3), increasing for k — 0.01, peaking when k = 0.05 and decreas- 
ing with k up to 0.5. At peak we found C = 3.11 ± 0.2, L = 0.62 ± 0.05 (Fig. [SjA.); these values were 
comparable to the average ones that we measure experimentally both on spiking and LFP activities and 
larger than those estimated on Erdos-Rcnyi networks {P < 0.001, ranksum test, Fig. [8]4_'). 
From our results we conclude that S has indeed a critical dependence on the number of functional 
networks. Overall, our experimental and simulation results suggest a sharp tuning in the number of func- 
tional networks that are active in cortical and subcortical nuclei. This number is likely to be preserved 
both in spontaneous and stimulated conditions. 



Discussion 



In this paper we show that spiking and LFP activities of neurons, distributed in contiguous or distant 
populations, present clear signs of small-world functional organization with sub-second invariance. Fur- 
thermore, we show that this distinctive functional organization persists in the presence of tactile stimuli, 
independently of their intensity. Our results further suggest the existence of efficient communication 
mechanisms, driven by LFP phase locking, between neuronal ensembles, distributed in distant brain ar- 
eas. These findings gain particular relevance in brain operations involving large co-activating neuronal 
populations engaging in demanding tasks 27 . In particular, in the thalamocortical loop our findings pos- 



tulate a solid and fast mechanism to arise synchronization response in sensory detection. Eventually, we 
confirmed these findings on simulated networks with hundreds of thousands of nodes (synthetic neurons) . 



Small-world in brain networks 



Studies in anatomical brain connectivity detected small-world networks 14 16 28 pointing out that 
this topology, well appreciable in humans, is shared also by animals with different degrees of brain 
development [9 12 29 . The complementary observations that brain pathologies like epilepsy 30 and 



schizophrenia 31 are well mirrored by signs of small-world network disruptions with abnormal functional 



connectivity further supports the interpretive and explanatory strength of small-word topology. 
Aside from these results obtained on macroscopic networks, few studies have focused on the functional 
connectivity at the neuronal level. The few studies available showed intrinsic small-world functional 
connectivity of single brain areas in the visual areas of primates 14 15 . These beautiful approaches 



give access to the richness of the local functional connectivity scaled at the neuronal level but do not 
provide, however, a description of connectivity framed on a wider extent and involving different neuronal 
populations nested in close or far brain regions. Just because of the composite organization of the 
somatosensory messages, the analyses of functional connectivity limited to the intrinsic networks of one 
single area would miss the rich repertory of multiple functional distributed concordances and contributing 
to their construction. From our results it seems that the functional organization in small-word networks 
is ubiquitous at different spatial and temporal scales. This organization represents a stable scheme for 
information processing in brain areas, networks and neurons. There are further features inherent to the 
dynamics of small- worldness. In addition, one of the most important properties is the number of nodes, 
namely the number of neurons involved in the interplay among the involved areas. The small-worldness 
can be expressed with a value that depends, for the small-worldness to exist in a network, on the number 
of the composing nodes. The statistics for capturing this essential trait (the S statistics, see Methods) 
shows that, for very large node populations, with thousands of nodes, as from in fMR imaging, the value 
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of S may be greater than 100 10 



For a minor node number (24 to 83 in our single unit recordings) , the 
in accordance with related literature 
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This 



value of small- worldness S lies between 1.2 and 4.5 
discrepancy between large and small networks and their relative small-worldness appears to be originated 
even in the node numerical disproportion. This hypothesis has been empirically proven by Barrat et al. 
who showed that small-world networks can exhibit large values of small-worldness only for sufficiently 
large network size [s]. 



Neuronal assemblies and LFPs 

Our findings support the definition of Hebbian neuronal assemblies as functional units organizing neurons 
distributed in several areas. The marker of such a functional connectivity is the transiently coordinated 
activity (both spikes and LFPs) found in our recordings. Furthermore, another important consequence 
emerges from our results: Canolty and colleagues suggest that neuronal oscillations enable selective and 
dynamic control of distributed functional cell assemblies [5]. We could enrich such scenario speculating 
that LFP coherent activities may drive the integration of local computations occurred in these distributed 
cell assemblies. Thus, the small-world topology expressed by synchronizing and distributed LFP phases 
supports a hierarchical and efficient functional substrate for incorporating neuronal assemblies. In a 
more summarizing view, the distributed LFP-to-spike interactions show that neuronal assemblies are not 
homogeneous populations but, potentially, a collection of diversely operating small- worlds incorporating 
the examined LFPs in the different contexts. 



Computational Model 

To confirm the validity of the small-world hypothesis on large networks, we extended the data on a 
computational model made of 100000 neurons embedding 5000 small-world networks. Some important 
assumptions taken into consideration in our network implementing must be described more in detail. 
Namely, this computational model holds on three fundamental assumptions: i) neuronal populations 
compute information, while perform behaviors, following a small-world functional organization; ii) indi- 
vidual neurons can partake to several of these functional populations (point 4 in structural assumptions, 
see Methods); iii) these functional populations can work in parallel (point 3 in functional assumptions, 
see Methods). The first claim, the working- hypothesis of this paper, is partially supported by our re- 
sults and by numerous suggestions in the recent literature [9 12 . The second claim is supported by the 



recent findings on dendritic computations 32 - 36 stating that different dendritic branches can partici- 



pate in distinct neuronal computations. Eventually, the last point is based on recent observations about 



neuronal interactions [5j[20 21 suggesting that LFP activity, that influences spiking activity, incorpo- 
rates many distinct phase signals simultaneously as a multidimensional extension of the two dimensional 
communication-through-coherence hypothesis. Moreover, small-world topologies are associated with high 
global and local efficiency of parallel information processing [10, , consequently if we assume small- world 
network as model for network functional organization, in the same manner, we must assume the intrinsic 
parallelism permitted by small-worldness. Even from this result we conclude that neuronal assemblies 
are functionally organized as strong small-world networks and the recorded weaker small-worldness can 
be justified by the limitation of recording and by the cooperation of each neuronal assembly in multiple 
tasks. 



Conclusive notes 

In an evolutionary perspective, small-world topologies appear preferentially selected among network 
topologies under the natural constraints (efficiency and reliability) of brain expanding complexity in 



the mammal phylogenetic tree 11 . It effectively satisfies the necessity of internal input integration and 



grants for best responses to environmental requirements. Moreover, small-world networks seem coherent 
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with the recent advancements in the physiology of neuronal networks. More punctually, from a functional 
point of view, small-worldness appears to provide dynamical features, such as communication-through- 
coherencc l^^Oj^^i for fast information integration where even far located nodes efficiently participate 
to the integration process. 



Materials and Methods 



Ethical notes 

All the animals have been treated along the Italian and European Laws on animal treatment in Scientific 
Research (Italian Bioethical Committee, Law Decree on the Treatment of Animals in Research, 27 Jan 
1992, No. 116). The National Research Council, where the experiments have been performed, adheres 
to the International Committee on Laboratory Animal Science (ICLAS) on behalf of the United Nations 
Educational, Scientific and Cultural Organizations (UNESCO), the Council for International Organiza- 
tions of Medical Sciences (CIOMS) and the International Union of Biological Sciences (lUBS). As such, 
no protocol-specific approval was required. The approval of the Ministry of Health is classified as "Biella 
1, 3/2011" into the files of the Ethical Committee of the University of Milan. 



Animal preparation and stereotaxis 

Seven male albino rats (Sprague-Dawley, Charles River, Calco, LC, Italy, 300-400 g) were chosen in the 
set of 11 animals employed in the recordings. All the animals were maintained in 16/8 hours daylight- 
dark regimes, food and water ad libitum. The rats underwent preliminary barbiturate anesthesia for the 
surgical experimental preparation. The jugular vein and the trachea were cannulated to gain, respectively, 
a drug delivery access and the connection to the anesthesia-ventilation device. Before the placement of 
electrodes the rats were paralyzed by intravenous Gallamine thriethiodide (20 mg/kg/h) injection and 
connected to the respiratory device delivering (1 stroke/s) an Isoflurane (2.5% 0.4 to 0.8 1/min) and 
Oxygen (0.15-0.2 1/min) gaseous mixture [37[ . Curarization was maintained stable throughout the whole 
experiment by Gallamine refracted injections (0.1 ml of the original solution/h). The anesthesia levels 
were maintained into ranges avoiding any corneal or retraction reflex (in absence of curarization) with 
low intensity noxious mechanical stimuli applied on a posterior paw [37 . 

We chose three areas for the simultaneous neuronal recordings in the left brain: the thalamic median 
nuclei (in particular the centromedial (CM)), the thalamic ventro-postero-lateral nuclei (VPL) ,18(|38] 
and the primary somatosensory (SI) cortex. Fast tactile stimuli were delivered to the right posterior paw 
(see Fig. [TJ^") by a suitable electromechanical device. 

Two holes were drilled on the skull respectively. A 3 x 2 mm bone window for the access of the cortical 
matrix electrodes and a larger bone window (6x2 mm) allowing for the simultaneous insertion of 
two parallel electrode matrices directed to the thalamic nuclei. The cortical access was set around 



a reference at -1.5 mm AP and -2.5 mm ML on the left, 39 and the electrode matrix was driven 
around 450 to 800 micrometer deep by an electronically controlled microstepper (Narashige, Japan). The 
thalamic access was centered at the two focus points of -6 mm AP, -0.8 and -2.5 mm ML ^39]. The 
electrodes were inserted with a 25° slant and driven at least to 5500 /im in depth and then advanced by 
a second electronically driven microstepper (AB Transvertex, Stockholm) until responses were observed 
to peripheral test stimuli. The neuronal recordings were obtained with two types of matrices, a vertical 
array devoted to the cortical recordings and two planar matrices devoted to the thalamic recordings. The 
vertical array was a multitrode (Multitrode Type 1, Thomas RECORDING GmbH, Giessen, Germany) 
with 8 gold contacts (125 fim contact spacing) with an average impedance of 1.2 Mil. For planar matrices 
were 3x3 frames of tungsten or Pt-Ir electrodes, inter-tip distance 150-200 fim, tip impedance 0.5-1 Mil 
(FHC Inc., ME, USA). 
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Fast thalamic and cortical responses to light tactile stimuli in the plantar aspect of the right hind limb 
were used as anatomo-functional acceptance criterion for acquisition. 



Tactile Stimulation 

Controlled stimulation was delivered through a blunted cactus thorn on each of five sites of the rat right 
hind limb (Fig. The tip was mounted on the dust cap of a speaker and driven through an Arduino 
microcontroller board (available at http://www.arduino.cc). At the beginning of each stimulation epoch 
the tip was lightly placed over the skin. Fast 5 ms pressure pulses were applied following a semi-random 
sequence. Pulses occurred in couplets. The delay between the first pulses of each couplet was set at 500 
ms. Every second pulse of each couplet followed the first by a random delay extracted uniformly in the 
range 150-250 ms (see Fig. S4C). The stimuli semi-randomness was adopted to avoid habituation. 



Neurophysiological recordings and preliminary data analyses 

For signal amplification and data recordings a 40 channel Cheetah Data Acquisition Hardware was used 
(Neuralynx, MT, USA, sampling frequency 32 kHz). Electrophysiological signals were digitized and 
recorded with bandpasses at 6 kHz/300 Hz for spikes, 180 Hz/1 Hz for Local Field Potentials. The data 
stored were analyzed off-line both using Matlab and by locally developed software. We chose a total of 
404 neurons (66 ± 17 in each experiment) out from the set of the acquired signals. The neural firing rates 
had a mean of 34.1 Hz with standard deviation of 26.8 Hz. Sorted cells with average rates below 4 Hz 
and above 100 Hz were excluded from the analysis. Distribution of firing rates and inter-spike intervals 
is shown in the Supplementary Figures S4A-B. After the recordings the LFPs were downsampled to 0.5 



KHz. We used for filtering the same techniques described in 40 . After filtering and downsampling, the 
spike contamination of LFP signals was null avoiding further spike removal techniques 41 . The spikes 
were extracted and sorted by using the Wave-clus MATLAB toolbox ^42j. The timestamps of spike 
occurrences were represented by binary sequences where I's labeled a spike. We considered time bins of 
1 ms thus avoiding occurrence of multiple spikes within the same bin. Finally, we split each sequence 
into fixed-length (from 50 to 1000 ms) overlapping windows (Fig. [ip), thus obtaining an ordered set of 
equal length sequences. 



Functional connectivity by spike-train similarities 

Interactions between neurons can be very complex, time-delayed and asymmetric. This can represent a 
problem in a experimental configurations where the communications between distant neurons are taken 
into consideration. Standard techniques like correlation analysis are, in many cases, unable to detect 
such events. Indeed dependencies between neurons from thalamus and cortex can last tens of millisec- 
onds ,43-45, . A toy example is shown in Figure [l]A_-B where an indirect (but actual) interaction exists, 
respectively, between neurons A and B (A ^ B). The neuron A is recorded and its activation triggers the 
long chain that, from site A, produces firing activity to neuron B in another brain region. The direct path 
between neurons C and B allows, instead, for synchronous spike patterns well detectable by correlation 
analysis (Fig. [lj3). 

In general, in simultaneous recordings from separated brain sites, it is unlikely to find couples of physically 
wired neurons although axonal projections connect the sites. It's definitely more probable that spiking 
activity relations may reflect a complex physical substrate, where chains of activations exist between 
them provoking the observed coherent activity. Such problem can be solved by mathematical tools able 
to model arbitrarily long relationships. 

In this work, we proposed a novel framework wherein spike trains with arbitrarily long temporal de- 
pendencies are modelled by Markov stochastic models. Typically, in regular Markov models, each state 
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depends only on the previous state while higher-order Markov models suffer from high state-space com- 



putational complexity 46 . Here, we used the Variable-order Markov Models (VMMs) 47 48 because 



they are able to address these limitations. Lossless compression algorithms (LCAs) represent one of the 



most efficient techniques to estimate VMMs 48 . Within the set of LCAs we chose the Prediction by Par- 
tial Matching (PPM) algorithm fio'.'SO^ which is considered the best match between prediction accuracy 
and speed [48j. In the last 15 years, Vitanyi and colleagues have developed a function, the Normalized 



Compression Distance (NCD) [51 -54 , which estimates the distance between symbolic sequences. This 
function performs the estimation directly by the sizes of the compressed sequences. In fact, the better is 
the VMM estimation the shorter is the compressed sequence size. 

In this work, we redefined the NCD function in order to reverse its assigned values pointing to similarity 
instead of distance. We called such function the Normalized Compression Similarity (NCS). Formally, 
given that x and y are two neural sequences (e.g. spike trains or LFP phases), the NCS is defined as 
follows: 

NCS{x, y) = l- NCD{x, y) = ^^^{c {x) , C {y)) 

where the C function represents the compressed sequence length and • is the sequence concatenation 
operator (e.g. 0101 • 101 — 0101101). If NCS{x,y) is close to 1, the sequences x and y are considered 
similar. If close to 0, the sequences are strongly dissimilar. We therefore evaluated the NCS function 
on time windows (50-1000 ms) of the recorded (binary) spiking activity (Fig. [ip-D) assuming that high 
values of similarity (greater or equal to 0.85) corresponded to actual functional connections (Fig. [ijil). To 
note that significant NCS values do not imply significant correlations but the opposite is true. Although 
the NCS is asymmetric function, it is not relevant for the causal interaction analysis and consequently 
for effective connectivity. 

Because NCS has never been used as method to estimate the functional connectivity in neurophysiological 



data (although many LCAs has been used as information measure of synaptic efficacy 55 and spike 



train similarities 56 57] ), our aim was to prove that NCS effectively captures similar asynchronous 
spike patterns in comparison to correlation coefficient. Furthermore, we demanded that NCS did not 
bias its estimations with respect to independent spike trains generated by simulations. To address these 
requirements we performed two experiments: i) we evaluated the NCS and the Pearson correlation 
coefficient over a set of independent binary spike trains and ii) we evaluated the NC S and the correlation 
coefficient with couples of binary spike trains (1000 ms long) containing an identical short and random 
spike pattern (50 ms) that is fixed and centered in the first train and drifts (from left to right) in the 
second trains (see Fig. S5B, overall 37 drifts). In the first experiment we found that the NCS method 
did not show any significance for independent spike trains as well as the Pearson coefficient (mean = 
-1.6 • 10-^ std = 0.005 for Pearson coefficient and mean = 0.013, std = 0.004 for NCS, Fig. S5A). 
This guarantees that the NCS function is unbiased to independent distributed binary spike trains. In 
the second experiment, the NCS proved its efficacy to detect long-range interactions where Pearson 
coefficient fails. In fact, during the whole shifting epochs, significance held while with Pearson coefficient 
appeared only when the two spike patterns were very spatially close (Fig. S5C, drift number 18-23). 

Functional connectivity by LFP phase synchrony 

LFPs are low frequency signals reflecting a wide range of synaptic events. In this work we investigated 
the synchrony of LFP phases originated in different recording sites during spontaneous and tactile evoked 
activities. We measured phase synchronies between two recorded LFP sequences (a; and y) by the following 
function 

^(^X,y) = (^eJ(^^9(H{x))-arg(H{yy))'^ ^2) 

where e is Napier's constant, H is Hilbert Transform, arg is the argument function and i is the imaginary 
unit. The Hilbert transform and the argument were computed with, respectively, the hilbert and the 



10 



cLngle Matlab functions 24 41 . When j{x,y) is equal to 1 (0), then x and y are perfectly synchronous 
(asynchronous) . 

Complex Brain Network 

By using the NCS and 7 functions, we estimated the functional connectivity of the recorded neuronal 
networks. We first split each recorded sequence into equal-length time windows (Fig. [T]) and then we 
computed the adjacency matrix for all neurons. The resulting matrices exhibited values in the unitary 
interval. We repeated the analyses with different window sizes from 50 ms to 1 s (fixing the maximum 
dependency order to half of the window size). A fixed threshold (typically set to 0.85 in all experiments) 
selected the strongest connections, thus allowing for the construction of the functional connectivity graphs. 
For instance, if x and y represent the spike sequences of two distinct neurons with NCS{x, y) = 0.92, then 
a connection x ^ y is established. During the analyses, other values (from 0.75 to 0.95) for threshold 
are used without significant changings in results. 

The functional connectivity extracted from extracellular recordings can be represented by graphs. For 
the analysis of these graphs, we introduced a set of common statistics from the Complex Network Theory 
able to detect possible matches between the extracted graphs and prominent topologies like small-world 
networks. A small-world network is generally obtained by evolving a basic ring lattice graph, where each 
node is connected to their K neighbors. The chosen neighborhood involves typically much less nodes 
than the total node number N {N ^ K ^ In(A^) ^ 1). The graph evolution is achieved by randomly 
adding and removing edges from the starting graph [6||7|. The resulting graph has many, typically 
small, quasi-complete subgraphs (cliques) where each node is connected to every other node in the clique. 
Furthermore, small-world networks exhibit short average distances between nodes. 

From a functional perspective, small-world networks can express two important information processing 
features: information integration and segregation [12[|59| . Functional segregation recruits specialized 
processing within densely interconnected nodes (cliques). Functional integration combines information 
processed in distributed nodes or cliques. These network properties can be measured by two statistics: 
the transitivity (C) and the characteristic path length (L) [6O]. The former measures how close the 
neighbors of a node are to being a clique. The latter estimates the average shortest path length in the 
graph, i.e. how much the nodes are accessible. 

Formally, let G be a graph represented by the adjacency matrix A = aij where i, j N = {1, . . . ,n} the 
set of nodes. We define the measure C as 



where ti is defined by 



\ X! °-ii°-'.hajh (4) 



2 

j,heN 



and ki is the degree of the node i. We also define the measure L as 

L^-Y ^^^^'^^'^'^'^ (5) 
n ^-^ n — 1 

where dij is the shortest path length between nodes i and j. Finally, it is possible to combine both 
measures in order to obtain a global small- worldness measure [S) as 
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where C" and are computed on random networks with the same connection density of the experimental 
ones. Graphs are considered small-world networks whether 5^1 and L/U" > 1. Finally we used a 
further network statistics called betweenness centrality, a function of a node that estimates the number 
of shortest paths from all vertices to all others that pass through that node. It is an indirect measure of 



the load of a node within the network 60 



Computational Model 

The validity of the results obtained from experiments has been further tested by computer simulations. 
Because of several limitations in in vivo experiments, we developed a computational model to verify 
observed evidences on large-scale networks. Furthemore, the model proposes a set of recent neurophysi- 
ological facts summarized in the following structural and functional assumptions: 
Structural assumptions: 

1. Each neuron is represented by a node. 

2. Functional connections between neurons are represented by edges. 

3. Neurons are functionally organized as small-world networks. 

4. Many small-world networks are embedded within the simulated network thus each node may belong 
to multiple small- world topologies (Fig. S6A-D). 

5. According to the typical amount of neocortical neurons in a microcolumn [26| , the sizes of small- 
world networks are randomly sampled by a uniform distribution ^ [/(50, iV/100) where N is the 
number of nodes. 

Functional assumptions: 

1. Each small- world network represents the processing of specific information (working hypothesis). 

2. At most two randomly chosen small- world networks are allowed to run independently in a time 
window. A single neuron can work either alternately or concurrently within different networks of 
pre and post-synaptic nodes. This assumption claims that neurons are not exclusive and that they 



can partake simultaneously to many (2 in our case) tasks 21 33 -36 



Neurons express their spikes within 10 ms time windows. Neurons fire following a centrality criterion 



based on the node betweenness centrality 60 . A neuron with high betweenness fires its spikes after 
a lower betweenness neuron. This dynamic picture respects the information integration-segregation 
paradigm of a hierarchical structure. Namely, neurons at lower layers fire before neurons at higher 
layers (Fig. S6E-H). 



4. 



Simulation timesteps are set to 1 millisecond. Spike propagation times are uniformly distributed 
within the range [1, 3]. 
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The algorithm govering this network is the follow: 
Input: n^swn, njnodes, timestep; 

Output: the binary runodes-hy-timesteps matrix trains 

0.05; 
fc 5; 

for i -s— 1 to n_swn do 

n -S— suhsaiwple (njnodes) ; 
G -S— random_watts_strogatz(7i, ; 
wins ^ computeWindows (G) ; 
foreach time windows in wins do 
if 0.01 < randO then 



centrality ^ BetweennessCentrality (G) ; 
foreach node A in centrality do 

trains[A, wins + centrality (A)] ^ 1; 

foreach output node B of A do 

trains[A, wins + centrality{A) + centrality {B)] -S— 1; 



end 
end 
end 
end 

where the function subsampleO selects a subset of nodes out of the set of available ones. The function 
rajidom_watts_strogatz returns a random graphs built following the Watts-Strogatz algorithm, the 
function computeWindows () computes the length of the execution windows for the current graph and 
returns the list of such windows. The function rsLadO returns a uniformly generated random number 
between and 1. At last, the function BetweennessCentrality () computes and returns the centrality 
of each node. 

The model has been developed using the Python environment [6l]. For the g eneration of small-word 
networks we used the networkx package (available at http:/ / net workx . lanl . gov/ *) . Fig. S6I shows a raster 
plot obtained by sampling the simulated activity. The darkest blue represents the no-spike event and 
the other spike colors are associated with the specific small-word topologies that generated them. The 
algorithm describes a core loop where each small-world network is first randomly created by the library 
routine watts_strogatz_graph (probability of rewiring equal to 0.05). Small- world networks take only a 
fraction of the total node number and several small- world networks can share a subset of their nodes. In a 
second stage, the generated small-world network expresses its spiking activity following the betweenness 
centrality as in point 2 of the dynamical assumption. Low level uniformly distributed noise is added 
to the spike propagation time. The spikes of the current network occur randomly in equal size time 
windows (10 timesteps). The spike activity is saved and the loop restarts (source codes are available at 
[http://code.google.eom/p/swn-neuronal-networks/ ). 




end 
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Figure 1. The proposed framework for the estimation of neuronal functional connectivity. 

(A) A recording session from thalamic and cortical regions. Arrows indicate the effective influence 
among neurons. The electrode tips record the neurons in dark red. (B) The firing patterns of the 
cortical neuron B produce common firing patterns both with neurons A and C but with different time 
delays. In particular, C — > -B (red spikes) can be easily inferred by correlation analysis instead of 
A B (blue spikes) hardly detectable. (C) Recorded signals are processed in overlapping windows 
lasting hundreds of milliseconds. (D) Spike trains are modeled by VMMs and compressed by LCAs. 
The functional connectivity strength between the spike trains A and B is estimated by the length of the 
compressed spike trains (C(A), C(B)) used by the NCS function. Whether NCS{A, B) is greater than 
a fixed threshold then we can conclude that A ^ B. (E) An example of functional graph extracted by 
recordings in the i*'' time window. (F) Typical sites of the rat paw for the tactile stimulation. 
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Figure 2. Small-world network statistics for spike functional graphs. The four Box- Whisker 
plots represent the computed statistics for graph with time windows of length respectively of 50, 250, 
500 and 1000 ms. Small-world properties are satisfied since S" > 1 in each plot. 
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Figure 3. Example of stimulus evoked redistribution of spiking functional connections. 

Functional weights are redistributed from the pre- (A) to post-stimulus (B) configurations. Red, blue 
and green nodes indicate neurons respectively from VPL, SI and CM. Some neurons, not functionally 
connected in the pre-stimulus graph, are involved in the stimulus processing (3-6,10,24,28,32,34,37,42). 
Conversely, some neurons, previously employed, are excluded in the functional graph (11). Furthermore, 
many neurons change their functional roles. For instance, cortical neuron 32 is a hub node (node with 
high degree) with many functional connections with VPL and CM thalamic neurons in (B) while is not 
involved in (A). Again, cortical neuron 23 constitutes a small clique with cortical neurons 22 and 40 in 
(A) while in (B) becomes a small hub node with diverse thalamic and cortical neurons. 
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Figure 4. Functional connectivity at spike level during stimulus evoked activity. A-D) The 
box-whisker plots sides represent the network statistics in pre- and post-stimulus conditions. (E-F) The 
Small-world network statistics remain stable over the five stimulation sites. 
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Figure 5. Small-world network statistics for LFP functional graphs. (A-D) The four 
box-whisker plots represent the computed statistics for graph with time windows of length respectively 
of 50, 250, 500 and 1000 ms. Difference matrices values for LFP phase. After stimuli the standard 
deviation of LFP evoked synchrony is significantly smaller in comparison with control values obtained 
in random sampled time windows. 




Figure 6. Example of LFP phases couplings. Functional connections are disposed on the 
LFP recording sites. Green nodes represent CM channels, red nodes represent VPL channels and 
blue nodes represent SI channels. (A) Before a tactile stimulation, LFPs are tightly coupled among the 
LFP channel of the same brain areas. (B) After an effective tactile stimulation, LFPs broke their 
inter-site synchronies and established cross-site phase couplings. 
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Figure 7. Small-world statistics for LFP phase sequences during evoked activity. (A-D) The 
box-whisker plots represent the network statistics in pre- and post-stimulus time windows. (E-F) The 
Small-world network statistics remain stable over the five stimulation sites. 
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Figure 8. Small-world statistics in the simulated model. Estimated statistics assuming multiple 
co-actives functional (A) small-world and (A) random topologies. (B) Small-worldness in comparison to 
the number of embedded small-world topologies. The ratio p indicates the level of embeddedness for 
computational models constituted by fixed number of nodes (100000) and variable embedded 
small- world networks (1000 to 50000). Small-worldness emerges for a particular range of admitted 
topologies (10 to 25). 
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Figure Legends 
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